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SUMMARY 



? 


This paper presents a numerical scheme, based on the finite element method, to solve strongly coupled 
fluid flow and heat transfer problems. The surface radiation effect for gray, diffuse and isothermal surfaces is 
considered. A procedure for obtaining the view factors between the radiating surfaces is discussed. The 
overall solution strategy is verified by comparing the available results with those obtained using this 
approach. An analysis of a thermosyphon is undertaken and the effect of considering the surface radiation is 
clearly explained. 


INTRODUCTION 

There are many engineering applications in which a coupled analysis of fluid flow and heat transfer 
is desired. Among a large list of such examples, a few important ones are design of heat exchangers, cooling 
of electronic components, climate control and underhood analyses in automobiles, performance of industrial 
furnaces, heat transfer analysis in confined cavities, and, cooling and heating of buildings, etc. The fluid flow 
analysis generally requires solution of conservation equations of mass and momentum. Several numerical 
approaches are available (refs. 1 to 4) under a variety of boundary conditions. In heat transfer studies usually 
energy conservation involving all three modes (namely, conduction, convection and radiation) is expected. 
However, until recently, conduction and convection heat transfer modes were accurately accounted for while 
approximations were made for including the radiation analysis (ref. 5). The high nonlinearity involved in the 
basic theory precluded from obtaining analytical solutions and a use of ordinary numerical methods for 
practical problems. The availability of cheaper computer resources has caught the attention of researchers 
wanting to include accurate radiation analyses in their studies. This is reflected in a collection of papers 
included in (ref. 6) published recently. 

The aim of this paper is to present a numerical methodology for analyzing fluid flow and heat transfer 
problems (including all three modes). A brief account of numerical solution of Navier-Stokes and continuity 
equations using the finite element method is presented. The assumptions involving the heat transfer via 
radiation include non-participating fluids and gray, diffuse surfaces based on enclosure theory (ref. 8). 

Solution of strongly coupled (heat transfer and fluid flow) phenomenon with natural convection is 
demonstrated through a couple of examples. To benchmark the developed code, a comparison with the 
already reported results is made. This is followed by a discussion of results in an analysis involving a study of 
thermosyphon (ref. 9). a passive system used for cooling of electronic components. 
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GOVERNING EQUATIONS 


In this section, the basic equations associated with the fluid flow and heat transfer are discussed. 
Generally, it suffices to consider the conservation of mass, momentum, and energy in the given domain of 
interest. In the presence of surface radiation, additional equation representing the conservation of radiative 
energy must also be considered. The effect of radiative fluxes on the relevant surfaces must be reflected in the 
overall energy balance. In summary, the following equations must be solved to conserve mass, momentum, 
energy, and radiative energy: 

conservation of mass: 
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conservation of momentum: 
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conservation of energy: 
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For explanation of the symbols employed, refer to the section titled Nomenclature. It should be noted that 
Equations ( 1 ) through (3) are used for incompressible fluid flow with Boussinesq approximations invoked to 
model the natural convection phenomenon. 

conservation of surface radiative energy: 
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In deriving Equation (4), it is assumed that the surfaces are gray, diffuse and isothermal (ref. 8). The view 
factors, Fij, between surfaces i and j, appearing in Equation (4) must be computed when attempting the 
solution of this equation. In the next section, a discussion on view factor calculations is undertaken. 


176 


COMPUTATION OF VIEW FACTORS 


In order to compute q r ’s (in Equation (4)), view factors Fjj, between all radiating surfaces must be 
available. In this section, the physical meaning of viewfactor and its calculation will be discussed. For a 
better understanding, ij in the above equation can be replaced with I and 2. Thus, view factor. F i_2, between 
two arbitrary surfaces (see Figure 1), ‘ 1 ‘ and ‘2’ is defined as a fraction of diffuse radiant energy leaving 
surface ‘ 1 ’ that arrives at surface ‘2’. Mathematically, 


F 1-2 = 
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where A| and A 2 are the areas of surfaces 1 and 2, respectively, ri2 is the distance between the two elemental 
areas dAi and dA2, ©l is the angle between the position dependent normal vector n*and the line connecting 
dAi and dA2. Angle 02 is defined in a similar way. It must be noted that Cos 0| and Cos 02 must be positive 
in order for the surface dAi and dA2 to 'see’ each other. If either of the cosines has a negative value, the 
corresponding view factor, FdA|-dA? should be set to zero. Such cases, in which the inactive side of the 
radiating face acts as an obstructer, will be termed as 'self-obstruction’ cases. Also, view factor F 1 -2 should 
be set to zero, if a third surface obstructs the view between surfaces 1 and 2. 

In order to calculate view factors internally, the user must specify the radiation surfaces in terms of the 
finite element faces of a discretized domain. The user must also specify which of the two sides is a radiatively 
active side. These pieces of information can be supplied very easily via the already existing card in the NISA 
file of NISA/3D-FLUID. Each radiating face is taken as one radiation surface. View factors between the 
radiating surfaces are automatically generated by NISA/3D-FLUID taking into account self-obstruction and 
obstructions due to a third surface. 

As can be assessed from the preceding discussion, computing view factors can result in usage of 
excessive computer time. To economize this computation, different techniques are used depending on 
whether the geometry being analyzed is 2D, 3D or axisymmetric. For example, double area integration 
method (ref. 8) is employed in comparison with contour integration method (ref. 8) when a 3D geometry, with 
radiation surfaces, is being analyzed. No special directives are required when computing view factors for 
axisymmetric geometries. NISA/3D-FLUID internally generates a complete 3D model (with the axis of 
symmetry as the X-axis [NISA/3D-FLUID]) to calculate the required view factors. Furthermore, for 2D 
problems, a completely different approach, called Hottel’s crossed-string method (ref. 8) is employed for its 
computational efficiency and accuracy. Reference 8 provides more details for evaluating view factors for 
interested readers. 


FE FORMULATION & SOLUTION PROCEDURE 

The partial differential equations (Equations 1 through 3) and the radiative balance equation (Equation 
4) are to be solved simultaneously to account for the fluid flow and heat transfer analyses in a given domain 
with specified boundary conditions. The convective terms appearing in Equations (2) and (3). simultaneous 
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solution of Equations (3) and (4), and arbitrary geometries encountered in most practical problems would 
require numerical tools for obtaining solution to coupled Equations 1 through 4. The Galerkin method in 
conjunction with the finite element method (ref. 8) form the basis of discretizing Equations 1 through 3. The 
penalty approach (ref. 3) is employed to eliminate the pressure from Equation (2) making use of Equation (1). 
For further details, refer to (refs. 3 and 10). The discretized form of Equations (2) and (3) can be written in 
matrix form as follows 

[KJ {X} = {1} (6) 

where Kij is the "stiffness" matrix, consisting of contributions from acceleration, diffusion and pressure 
gradient terms of Equation (2) and acceleration and diffusion terms of Equation (3). Xj represent [U, V, W ] 
for momentum equations and [T] in the case of energy equation. The vector fj is discussed more at length as 
this contains coupling terms in Equations (2), (3), and (4). For example, the vector {f} for the momentum 
equations is 

-J N 1 pgiP (T-Tp) dQ + j N 1 ( - P Sjj + tjj) njdT (7) 

n r 

Equation (7) indicates the influence of temperature distribution on the momentum equations while convective 
terms (included in Kij for Equation (3)) represent a dependence of the temperature field on the velocity 
distribution. 

Furthermore fj for the energy equation consists of the following term: 

( 8 ) 

f;= f N 1 Q d Q + f N 1 q d T 
J J n T 


where 


q = qa + qc + qr 


(9) 


In the above equality, q a , qc, and q r refer to the applied heat flux, effect due to convection boundary 
conditions, and that due to radiation on the boundary , respectively. The gray-body radiative effects can be 
considered via q r which is evaluated using Equation (4) for a "known" temperature distribution. It is thus 
evident that Equations (3) and (4) are coupled via q r and T. Far a complete enclosure. Equation (4) can be 
represented in the matrix as 

( 10 ) 

[R] {qr} = [S] {T} 
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Ideally Equations (6) for momentum and energy equations together with the radiative balance Equation (10) 
must be solved simultaneously. For practical reasons (computer memory and time, and nonlinearity in 
Equations (6) and (10)), a sequential approach is undertaken to solve these algebraic equations. Depending on 
the nature of coupling (strong for flows with free convective effects and weak for flows with forced 
convective effects), momentum, energy, and radiative balance equations are solved. For more details, refer to 
(ref. 10). It has been observed that q r (and hence T) solution may not converge or may do so slowly. An 
under relaxation of q r leads to its stabilization. This is achieved as follows: 

q[. + ' = «qr +l + ( 1 -a) qL (13 ^ 


where a is a user-defined relaxation factor. During a calculation sequence convergence checks are performed 
for velocity, temperature and surface flux, q r , distributions by evaluating the L 2 norms. The sequential 
calculations are performed until the L 2 norms of all the nodal variables and surface radiation fluxes fall below 
a user-defined tolerance. 


Special Cases: 

There are a few special cases which require a slight modification to the above methodology for 
including the gray surface radiative effects in the heat transfer analysis. These are as follows: 

a) Domain with plane(s) of symmetry 

b) Exchange of radiative flux through "windows" in the domain 

c) Exchange of radiative flux between the domain and surroundings 

d) Radiative surfaces with no thickness. 

The details of these modifications are presented in (ref. 10). 
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ILLUSTRATIONS 


The aim of the present paper is to discuss an efficient solution strategy that must be undertaken to 
solve coupled fluid flow and heat transfer problems in the presence of radiative energy exchange between 
gray surfaces in a domain of interest with specified boundary conditions. In the previous sections, the 
pertinent differential equations and their respective discretized forms (using the finite element method) are 
discussed. In this section, a discussion of the results obtained with the outlined procedure for a couple of 
problems is undertaken. 


Example 1 : Natural Convection with and without Surface Radiative Effects in a Cavity. 

The validation of the developed procedure is established by solving a problem studied by Behnia et al. 
(ref. 1 1 ). The fluid flow due to natural convective effects in a square cavity with radiating surfaces is 
considered. Figure 2 shows this cavity of a characteristic dimension, L and the specified boundary conditions. 
The top and bottom walls are adiabatic. The left wall is maintained at a uniform hot temperature, Th. The 
right wall has convective and/or radiative boundary condition. The convective heat transfer coefficient is h. 
The temperature of the surroundings and the ambient temperature are taken to be T«. All the internal surfaces 
of the cavity have an emissivity of 0.9 and the fluid in the cavity is air. The cavity size, L, can be chosen to 
get a Rayleigh number of 3xl0 5 . Table 1 shows a summary of conditions under which each case is analyzed 
with an aim of obtaining steady state temperature and fluid flow distributions in the cavity. Due to the 
presence of natural convective effects, strong coupling between the fluid flow and temperature fields is 
expected. The cavity is discretized into a graded mesh of 44 x 36 linear quadrilateral elements. The steady 
state algorithm of the code is invoked. Table 2 shows the relaxation parameters employed for each of the run 
detailed in Table 1 and the corresponding numbers of iterations required to obtain converged solutions. 

Figure 3 shows the isotherms obtained for the cases denoted as R300, EC300, and REC300. A 
comparison of isotherms for these cases clearly indicates the effect of surface radiation on the adiabatic walls 
(top and bottom), the isotherms are no longer normal to these walls. Figure 4 shows the streamlines for the 
cases R300, EC300, and REC300 respectively. Table 3 shows a comparison of the maximum value of stream 
functions obtained for these runs with those listed in Behnia et al (ref. 1 1). A good quantitative agreement 
between the results is evident. Figure 5 shows the horizontal velocity along the vertical center line for these 
cases. The velocity profiles shown in the figure compare well with those in Figure 7 of ref. 11. 

Example 2: Analysis of a Planar Thermosyphon. 

In this example, the fluid flow and temperature distributions are studied in a thermosyphon including 
the surface radiative effects. A thermosyphon is a device used for cooling of electronic components, heat 
removal systems for nuclear reactors, and having applications in solar systems (ref. 9). Since thermosyphons 
involve no blowing or pumping of fluids, they are less expensive and more durable (termed as passive 
systems) as these do not require external signals for operation. A schematic of planar thermosyphon and the 
assigned boundary conditions is shown in Figure 6. An analysis of fluid flow and heat transfer in a 
thermosyphon is presented in (ref. 9) without the surface radiation effects. These effects have been included 
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in the study here. All the walls are black and are considered to be radiating. By observing the radiation 
surfaces in Figure 6, it is evident that all the surfaces cannot "see" each other. In other words, the view 
factor computation, in th$ presence of third surface obstructions, is invoked. These computations are more 
complex and handled efficiently in NISA/3D-FLUID (ref. 10). 

First, the results are presented for the case in which only the convective effects (due to natural 
convection) are considered. The same analysis is performed in (ref. 9), in which the effects of varying the 
Rayleigh no. and a ratio of thermal conductivities of solid to fluid are considered. Therefore for the sake of 
comparison, results are presented for a Rayleigh no. of 10 4 and a ratio of thermal conductivities of 1 (see (ref. 
9) for more details). Figure 7 shows the stream function distribution for this case and the corresponding 
isotherms are shown in Figure 8. A good agreement between these results and those presented in (ref. 9) is 
observed. 

Now, the surface radiation effects due to the surfaces shown in Figure 6 is considered. The results 
for this case are not presented in (ref. 9). Figures 9 and 10 show distributions of stream functions and 
isotherms. A comparison of isotherms shown in Figures 8 and 10 indicate a considerable difference in their 
distributions. A further comparison of the velocity distributions, Figure 1 1 , at "inlet" and "outlet" of the 
thermosyphon show marked differences. The difference in these velocity distributions amounts to a 
difference of 25% in flow rate. This analysis clearly indicates that if the surface radiation heat transfer is not 
accounted for, inaccurate distributions of temperatures and velocities may result. 


CONCLUSIONS 

A numerical scheme based on the finite element method is presented for solving coupled fluid flow 
and heat transfer problems in the presence of surface radiation. A sequential solution of momentum, energy 
and, radiative energy equations is considered for efficient computer memory management and disk usage. 

The computed results validated the numerical procedure adopted for an analysis of coupled fluid flow and 
heat transfer phenomena. The results presented compared well with those reported in literature. It is shown 
via the results discussed in this paper that the surface radiative effects must be considered for a complete heat 
transfer analysis. More research is underway to extend this work to consider non-gray surfaces and eventually 
participating fluids. 
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Figure 5 Horizontal Velocity Distributiuon Along the Vertical Centerline for Runs 
REC300, EC300 and RC300 



Figure 6 Schematic Diagram for Example 2 (Thermosyphon) 
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Table 1 Summary of Three Different Runs 


Run 

Ra 

B.C. on the right wall 

Surface Radiation J 

REC300 

300,000 

Convective + radiative 

included 

EC300 

300,000 

Convective 

not included 

R300 

300,000 

Radiative 

included I 


Ra = GfPr ■ 

M 2 k 


Table 2 Relaxation Parameters and Number of Iterations 


I Run 

Relaxation Parameter 

No. of iterations 

Velocity 

a 

Temperature 

Y 

Radiative heat 
flux 
O 

REC300 

0.04 

1.0 

0.1 

64 

EC300 

0.04 

1.0 

0.1 

40 | 

R300 

0.04 

1.0 

0.1 

68 1 


Table 3 Values of |4"| ma x for Different Runs 


Run 

,u,n I'PIwm k 

ITlmax 3 * , a = — - 

a pCp 

Behnia et al. 
(ref. 11) 

NISA/3D-FLUID 
(ref. 10) 

REC300 

13.04 

12.94 

EC300 

10.93 

11.01 

R300 

11.93 

.1.79 
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